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In this paper, we describe a numerical approach to evaluate Feynman loop integrals. In this 
approach the key technique is a combination of a numerical integration method and a numer- 
ical extrapolation method. Since the computation is carried out in a fully numerical way, our 
approach is applicable to one-, two- and multi-loop diagrams. Without any analytic treatment it 
can compute diagrams with not only real masses but also complex masses for the internal parti- 
cles. As concrete examples we present numerical results of a scalar one-loop box integral with 
complex masses and two-loop planar and non-planar box integrals with masses. We discuss the 
quality of our numerical computation by comparisons with other methods and also propose a 
self consistency check. 

1. Introduction 

A method for the systematic calculation of loop diagrams is required to get precise theoretical 
predictions for elementary particle interactions. In this paper we present a fully numerical 
approach toward this computation. For simplicity we consider scalar loop integrals throughout 
the paper. The general expression of scalar loop integrals in Feynman parametric representation 
is 

where L is the number of loops, N is the number of internal particles and n is the num- 
ber of space-time dimensions. Here, C and D are polynomials of the Feynman parameters 
(Xj, i = 1, - ■ ■ , N) and they are determined by the topology of the corresponding diagram. 
An infinitesimal parameter, e, is introduced to make the denominator non-zero throughout the 
integration domain. With n = 4, the form of the expressions of the integrand for each and L 
is shown in Tabled] 
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Table 1 

Integrand of scalar loop integrals for the case n = 4 up to L = 2 and N = 7. For L = 1, C = 1. 
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A standard (analytic) method for multi-loop integrals is by a reduction to a set of integrals 
using integration by parts However, reduction often yields a large number of loop integrals 
and it becomes difficult to obtain an accurate numerical result due to stability problems and 
cancellation error, with a loss of trailing digits. On the other hand, in our approach it is not 
necessary to reduce the integral and the number of integrals to be computed is limited. 

Generally speaking the numerical computation of loop integrals becomes harder with an in- 
creasing number of loops and/or external legs since the behavior of the singularities becomes 
more complicated. We use an automatic integration technique with an efficient numerical ex- 
trapolation method and, if necessary, a suitable variable transformation of the Feynman param- 
eters. Since all computation is completely numerical, it is trivial in our approach to extend the 
number of loops and legs, be it at the price of an increased amount of work as measured in 
the number of integrand evaluations. Another advantage of the general numerical approach is 
that it does not matter whether the masses of the internal particles are real or complex. So far 
we have applied our approach to the computation of one-loop vertex (L = 1, = 3), box 
(L = 1, A^ = 4) and pentagon (L = 1, A^ = 5) as well as two-loop self energy (L = 2, A^ = 5) 
and vertex (L = 2, A^ = 6) diagrams [|2l[3l|4l[51[6l|7l. Numerical results by our method show 
good agreement with that of other methods [I|9l[iailllll2[ll[13[13[ia[l7l[l8l[ia|20l|2ll. 

Subsequently in section [2] we give a brief explanation of our approach with respect to the 
numerical techniques. In section [3] we present results of one-loop and two-loop diagrams with 
real or complex masses up to L = 2, A^ = 7 as examples. 

2. Direct Computation Method 

As in Eq. ([T]), ie is introduced in the denominator of the integrand. For instance, for the case 
L = 1 and A^ = 3 in Table [B we separate the real and imaginary part of the integrand as 

D — le + 

'^m ^ , = / (3) 

since our numerical integration package Quadpack [|22l currently does not support complex 
integrands. While analytically e is thought of as infinitesimal, we will replace it by a sizable 
number in a sequence of the form e = = a^'-^^^j = 0, 1, ■ ■ ■ where a is a positive number 
and / is an integer. When ej is finite, the integral converges and numerical methods can be 
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applied to the integration of Eq. For instance, varying ej geometrically as l.2^^^~^^ with 
a = 1.2 and / = 30, we get a sequence of integrals /(cj) corresponding to each e^ . It is our goal 
to obtain the limit limej-+o li^j)- This is an extrapolation as ej — )■ and we can accelerate the 
convergence of the sequence by an appropriate acceleration technique. We refer to this method, 
i.e., the combination of the integration and the extrapolation, as a Direct Computation Method. 

2.1. Multi-dimensional integration 

Basically in our method, we can choose any numerical integration procedure if it gives the 
numerical result to enough accuracy. However, most numerical integration methods fail due to 
singularities in the integration domain. We use the DQAGE routine from the one-dimensional 
Quadpack package [ 22] for a repeated integration [|2l[3l|23'| in the coordinate directions of the 
multi-dimensional integral. DQAGE is an adaptive quadrature routine. Generally it will partition 
around an integrand singularity "hot-spot" within the integration region, for an arbitrary location 
of the singularity. On each subinterval generated in the subdivision, DQAGE applies a variant of 
Gaussian quadrature where the sampling points are chosen by a Gauss-Kronrod scheme. 

2.2. Extrapolation 

As for the extrapolation, we choose an extrapolation method which does not require explicit 
information to be supplied about the sequence. Throughout this paper we present results using 
Wynn's e algorithm [ 24] as the extrapolation method. This works very efficiently even for 
sequences of /(e^) with slow convergence (providing the progression of tj is geometric). 

3. Computation of the loop integrals 

Here we show some numerical results by the Direct Computation Method. Since our approach is 
based on the complete numerical technique, the computation is possible for loop integrals with 
arbitrary masses no matter whether they are real or complex. The first example, in section [311. 
is a scalar one-loop box diagram with complex masses. In section 13. 2[ we treat scalar two- 
loop planar and non-planar box diagrams with real masses. For these examples we not only 
present numerical results but also discuss a technique for a cross-check of our numerical results. 
There are several ways for cross-checking, and comparisons with results by other methods 
are effective. However, if no other results are available for comparison, we propose a self 
consistency check to test the quality of the computation. 

3.1. One-loop box with complex masses 

In this section we consider a scalar integral /(s, t) for the one-loop box diagram shown in 
Fig.OC^ = 1, iV = 4 in Tabled defined as 



I{s,t)= dx dy dz— — , (4) 

Jo Jo Jo [D-iey 

where D is given by 

D = pIx^ + ply'^ + tz^ + {p\ +pI- s)xy + {pj -pl + t)xz + {pi -pl + t)yz 
+ + "^1 ~ rn^x + {—p\ + m\ — m^)y + {ml — m\ — t)z + m\, 

with s = {P1+P2Y = {P'a+PaY and t= {pi+Pa)'^ = {P2+PiY- One-loop box integrals have 
been fully analyzed and several powerful tools such as FF [|25l and LoopTools [|26l| have 
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been developed for their numerical evaluation. However, it is often tedious to include complex 
masses for the internal particles of the box diagram analytically. 
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Figure 1 . One-loop box diagram 



In 2007 and 2008 we compared our results for the one-loop box diagram contribution of gg — )■ 
hhH, to those obtained analytically by L. D. Ninh et al. {2JI- A severe numerical instability due 
to a Landau singularity was reported by L. D. Ninh for 211 GeV > Mr > 2Mw and 457 GeV > 
-\/s > mt in his numerical evaluation; and we observed the same instability using the Direct 
Computation Method. To regularize the singularity we, as well as L. D. Ninh and co-authors, 
included two complex masses for the internal particles as — initVt and — iMwYw, with 
Tt = 1.5 GeV and Tw = 2.1 GeV, respectively. After inclusion of the widths the instability 
disappeared and the results by L. D. Ninh showed very good agreement with ours [|7|. 

Subsequently in 2010 we computed the one-loop box integral with complex masses set to 
mj = 20 - Oi, ml = 10 - 5i, m| = 40 - lOi and mj = 10, 100, 1000 with pj = -60, 
pI = 10, pI = —10, pI = —10, s = 200 and t = —10. The results for both the real part and 
the imaginary part are compared with those by XLOOPS-GiNaC [ 28| by H. S. Do and P. H. 
Khiem. The detailed explanation was presented by H. S. Do in CPP 2010 [[30l|, showing not 
only their results and ours, but also results byLoopTools2.5 with the code DOC developed 
by D. T. Nhung and L. D. Ninh et al. [|29l|. All the results are in good agreement. 

Through the reported experiences we validated the Direct Computation Method for these loop 
integrals with real or complex masses of the internal particles. 

3.2. Two-loop box integrals with masses 

Here we discuss the computation of the two-loop box diagram for L = 2 and = 7 in Table[TJ 
Corresponding diagrams are shown in Fig. [2] and in Fig.[3l In the following let us consider the 
scalar loop integral defined as 



I(s,t) = — dxi dx2 dxs dx4 dx^ dx^ dx-j 5{1 — > xA-— (5) 
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Figure 2. Two-loop planar box diagram Figure 3. Two-loop non-planar box diagram 

D and C in the integrand for the two-loop planar box in Eq. ([5]) are given by 

7 

D = C x^m^ 

- {s{xiX2{Xa + Xs + Xg + X7) + X5X6(Xi + X2 + X3 + X4) + X1X4X6 + X2X4X5) 
+ tX3X4X7 

+ p\{x'i{xiXi + X1X5 + xixe + X1X7 + X4X5)) 
+ p\{x':i{x2Xi + X2X5 + X2X6 + X2X7 + X4X6)) 

+ p'i{x7{xiXi + X1X5 + X2X5 + X3X5 + X4X5)) 
+ pl^X-j^XiXii + X2X4 + X2X6 + X3X6 + X4X6))}, 

and 

C = (xi + X2 + X3 + X4)(x4 + X5 + xe + X7) - X4. 
Furthermore, D and C in the integrand for the two-loop non-planar box are given by 

7 

D = C Xf 

— {s(XiX2X4 + X1X2X5 + X1X2X6 + X1X2X7 + X1X5X6 + X2X4X7 — X3X4X6) 
+ t{Xs{-X4X6 + X5X7)) 

+ Pi{x3{XiX4 + X1X5 + XiXq + X1X7 + X4X6 + X4X7)) 
+ p\{xz{x2Xi + X2X5 + X2X6 + X2X7 + X4X6 + X5X6)) 

+ P3(xiX4X5 + X1X5X7 + X2X4X5 + X2X4X6 + X3X4X5 + X3X4X6 + X4X5X6 + X4X5X7) 
+ p1{xiX4Xq + XiXqXj + X2X5X7 + X2X6X7 + X3X4X6 + X3X6X7 + X4X6X7 + X5X6X7)}, 

and 

C = (Xi + X2 + X3 + X4 + X5)(xi + X2 + X3 + X6 + X7) - (xi + X2 + X3)^. 

For both diagrams we have that s = {pi + ^2)^ = (Ps + Pa)"^, ^ = (Pi + P^Y = {P2 + PaY and 
Pi + P2 + P3 + P4 = 0. 
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Since the number of dimensions of the integration is 6 and the behavior of the integrand 
is very complex, we performed a variable transformation of the Feynman parameters before 
numerical integration. With a suitable transformation for each diagram [ 33], we completed the 
computation of the real part of the two-loop planar box integral with mi = m2 = = uiq = 
m and iris = 7724 = rrij = M; and the real and imaginary part of the non-planar box integral with 
mi = 7712 = = ^6 = ITT- ^nd m3 = ms = rrii = M. For both pi = P2 = P3 = pI = 
The numerical results for each diagram are plotted as a function of fs = s/m? in Fig. |4] and 
in Fig. [51 respectively. For both cases, the kinematical parameters are t = — lOO.O^GeV^, 
m = 50.0 GeV and M = 90.0 GeV. 
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Figure 4. Two-loop planar box: Numerical results of real part in units of 10^^^ GeV^® for 
5 < /s < 25. Marks (red *) present results by the Direct Computation Method and marks 
(green □) are results by the reduction method. 



We cross-checked the numerical computation in two ways. First we made a comparison with 
the reduction method [i31J. After reduction we get numerical results using the Monte Carlo 
integration package BASES [fS^j. For the two-loop planar box both the results by the Direct 
Computation Method and those by the reduction method are plotted in Fig. HI The agreement is 
very good in the range of interest, 5 < < 25. 

For the two-loop non-planar box, both the real part for the whole range and the imaginary 
part for the range 4 < < 10 are plotted in Fig. |51 For the cross-check we attempted the 
reduction method. However, the behavior of the convergence by Monte Carlo integration is not 
very good for the intended range. Therefore six numerical results by the reduction method are 
plotted with error-bars in Fig. |5l Instead of the reduction for the cross-check, we performed 
a self consistency check using the dispersion relation between the real part and the imaginary 
part, 

= (6) 

■n J s — s' — le 
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Figure 5. Two-loop non-planar box: Numerical results of real and imaginary part in units of 
10~^^ GeV~^ for — 6 < < 10. Marks (red *) are results for the real part and marks (green 
□) for the imaginary part. Marks (light blue □) with error-bars present results by the reduction 
method. Marks (dark blue □) give numerical results by the dispersion relation for 5 < /s < 10. 



The results for the real part constructed numerically from the imaginary part using Eq. ^ are 
also plotted in Fig. [5] for the range 5 < < 10. 

4. Summary 

In this paper, we presented a complete numerical approach for multi-loop integrals. From a 
technical point of view, it is based on a combination of numerical integration and numerical 
extrapolation and we call the method a Direct Computation Method. Some numerical results 
are given as examples of the technique. Furthermore, for each numerical result a cross-check 
has been presented. We demonstrated the applicability of our approach to loop integrals of the 
one-loop box diagram with both real and complex masses and two-loop planar and non-planar 
box diagrams with masses. 
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